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ABSTRACT 

The  admissible  region  is  defined  as  the  set  of  physically  acceptable  orbits  (i.e.,  orbits  with  negative  energies).  Given 
additional  constraints  on  orbital  semi-major  axis,  eccentricity,  etc,  the  admissible  region  is  further  constrained,  result¬ 
ing  in  the  constrained  admissible  region  (CAR).  Based  on  known  statistics  of  the  measurement  process,  in  this  paper 
we  replace  hard  constraints  with  a  probabilistic  representation  of  the  admissible  region.  This  results  in  the  probabilis¬ 
tic  admissible  region  (PAR)  that  can  be  used  for  orbit  initiation  in  Bayesian  tracking.  While  this  is  a  general  concept 
that  is  applicable  to  any  measurement  scenario,  we  will  illustrate  the  idea  using  a  short-arc,  angles-only  observation 
scenario. 


1.  INTRODUCTION 

As  new  optical  sensors  come  online  and  more  and  more  optical  observations  become  available  for  space  objects 
previously  too  small  or  too  far  away  to  detect,  the  space  surveillance  community  is  presented  with  the  computationally 
challenging  problem  of  generating  initial  orbit  solutions  for  a  large  number  of  short-arc  line-of-sight  observations.  In 
order  to  perform  any  kind  of  probability -based  analysis  with  these  orbit  solutions,  we  require  an  accurate  representation 
of  their  uncertainty.  Properly  characterizing  the  uncertainty  will  allow  us  to  more  efficiently  deal  with  large  sets 
of  sparse  data  by  enabling  the  use  of  rigorous  probabilistic  techniques  to,  for  example,  perform  data  association, 
determine  collision  probabilities,  or  initialize  a  Bayesian  estimation  scheme.  This  paper  deals  with  the  problem  of 
characterizing  uncertainty  in  a  constrained  admissible  region  (CAR)  approach  to  initial  orbit  determination  (IOD). 

The  admissible  region  approach  has  been  previously  used  in  the  asteroid  tracking  community  by  Milani  et  al.  [1,2] 
to  deal  with  the  problem  of  identifying  asteroids  based  on  very  short  arc  observations.  Specifically,  they  referred  to  a 
region  in  the  plane  of  possible  ranges  and  range-rates  defining  those  values  for  which  a  given  line-of-sight  observation 
produces  an  orbit  solution  that  satisfies  certain  criteria.  This  concept  has  been  extended  in  the  SSA  community  by 
authors  including  Tommei  [3,4],  Milani  [5],  Farnocchia  [6],  Rossi,  Maruskin  [7],  Scheeres,  Alfriend,  Fujimoto  [8,9], 
Gronchi  [10],  Dimare,  Schildknecht,  Jehn,  DeMars  [11, 12],  Jah,  Schumacher,  and  Siminski  [13, 14]  to  deal  with  the 
problem  of  tracking  space  objects  in  Earth  orbit,  for  which  the  CAR  refers  to  a  region  in  the  range,  range-rate  plane 
which  produces  orbit  solutions  with  orbit  elements  satisfying  some  specified  bounds.  In  previous  work,  Schumacher, 
Wilkins,  and  Roscoe  [15, 16]  extended  this  concept  to  include  regions  in  the  range,  range  plane  satisfying  orbit  element 
bounds  for  pairs  of  observations.  Hussein  et  al.  [17]  applied  probabilistic  techniques  to  determining  the  admissibility 
of  an  uncertain  candidate  orbit  in  the  CAR  using  both  Unscented  Transform  (UT)  and  Monte  Carlo  (MC)  methods. 
Worthy  and  Holzinger  [18]  incorporated  measurement  uncertainty  into  the  admissible  region  approach  for  uncorrelated 
detections. 

An  unresolved  topic  in  the  SSA  literature  is  how  to  best  sample  the  CAR  to  perform  data  association  or  track  initiation. 
For  example,  Tommei,  Milani,  and  Rossi  [3]  use  a  Delaunay  triangulation  approach,  while  DeMars  and  Jah  [12]  (along 
with  a  number  of  other  authors)  use  a  uniform  distribution  approximated  as  a  Gaussian  Mixture  Model  (GMM). 
Simiski  et  al.  [13, 14]  use  an  Iso-Energy  grid  method  to  discretize  the  CAR  along  lines  of  equal  energy  (i.e.,  semi¬ 
major  axis).  The  goal  of  the  present  work  is  to  define  a  probabilistic  form  for  the  CAR  to  be  used  for  initializing  a 
Bayesian  tracking  scheme.  Based  on  known  statistics  of  the  measurement  process  and  assumed  statistics  in  the  semi¬ 
major  axis  and  eccentricity,  the  probabilistic  admissibility  region  (PAR)  can  be  determined  using  a  MC  method.  This 
resulting  PAR  which  is  made  of  a  set  of  particles  can  then  be  converted  to  a  GMM  using  the  Expectation-Maximization 
(EM)  algorithm  [19,20]. 
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The  paper  is  organized  as  follows.  In  Section  2.,  we  review  how  the  deterministic  CAR  is  constructed.  In  Section 
3.,  we  describe  how  to  construct  the  PAR  from  known  measurement  and  RSO  population  statistics.  In  Section  4.,  we 
describe  how  to  derive  an  MC  representation  of  the  PAR  and  in  Section  5.  we  describe  how  to  convert  this  particle 
cloud  into  a  GMM  using  a  versatile  and  robust  EM  algorithm  called  the  FJ-EM  algorithm.  We  provide  a  numerical 
example  in  Section  6.  and  conclude  the  paper  in  Section  7. 


2.  THE  ANGLES-ONLY  DETERMINISTIC  CONSTRAINED  ADMISSIBLE  REGION 

In  this  section  we  introduce  the  deterministic  CAR  for  short-arc,  angles-only  optical  observations.  In  the  next  section 
we  introduce  the  PAR.  The  starting  point  is  a  series  of  n  short-arc  right  ascension  and  declination  observations  («,.  6, ) 
taken  at  measurement  times  fj,  i  =  0, . . . ,  n  —  1,  with  epoch*  at  fo-  Following  the  second-order  polynomial-in-time 
least  squares  procedure  described  by  Maruskin,  Scheeres,  and  Alfriend  [7],  one  can  obtain  an  estimate  of  the  optical 

attributable  vector  (do,  do,  c>o,  c>o)  at  fo  (we  omit  details  here  for  brevity).  For  the  rest  of  this  paper  we  will  drop  the 
tilde  and  subscripts  for  ease  of  notation. 

We  now  review  how  to  construct  the  CAR  from  these  estimates  of  the  angles  and  their  rates.  We  mostly  follow  the 
notation  and  derivation  given  by  DeMars  and  Jah  [12]  (with  minor  corrections).  Let  r  be  the  inertial  position  of  the 
object,  q  be  the  inertial  position  of  the  ground  station  and  p  be  the  inertial  position  of  the  object  with  respect  to  the 
station.  Hence,  we  have  following  position  and  velocity  relationships 

r  =  q  +  p  (1) 

and 

r  =  q  +  P-  (2) 

Using  spherical  coordinate  system  with  respect  to  the  ground  station,  p  and  its  inertial  time  derivative  can  be  expressed 
as 


P  =  PUp 

p  =  pup  +  pa  cos<5uQ  +  pSus,  (3) 

where  up,  uQ  and  u,,  are  the  spherical  coordinate  system’s  basis  unit  vectors  and  are  given  by 

up  =  (cos  a  cos  S,  sin  a  cos  S,  sin  5) 
u„  =  (-  sin  a ,  cos  ct,  0) 

=  (—cos  a  sin  <5,  —  sin  a  sin  6,  cos  S) . 


The  position  and  velocity  r  and  r  are  related  to  the  two-body  energy  via  the  equation 

ll  ||2 


After  expressing  r  and  r  in  terms  of  the  parameters  ( a ,  d,  5,  S,  p.  p)  this  equation,  it  can  be  shown,  is  given  by 

p 2  +  w\p  +  F{p)  —  2£  =  0 


(4) 


(5) 


where 


F{p)  =  W2p2  +  W3p  +  W4 


2/i 

a/p2  +  w5p  +  w0 


*  Other  epochs  can  be  considered  as  well.  For  example,  the  mid-point  of  the  arc  may  reduce  error. 
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and  where  the  parameters  Wi  (i  =  0, . . . ,  5)  are  given  by 


w0  =  ||q||2 
wi  =  2q  •  up 
u>2  =  d 2  cos2  8  + 

W3  =  2 d  cos  (5(q  •  uQ)  +  2  j(q  ■  ud-) 
w4  =  ||q||2 
w5  =  2q  •  up 


Given  a  value  for  £,  one  can  solve  for  the  corresponding  curve  in  p—p  space.  The  curve  corresponding  to  £  =  0  defines 
the  admissible  region  inside  which  all  orbits  with  negative  energy  (£  <  0)  exist.  It  is  often  the  case,  furthermore,  that 
we  have  additional  inequality  constraints  that  one  can  impose  in  order  to  further  limit  the  set  of  possible  orbits  with 
negative  energy.  Two  important  such  constraints  are  constraints  on  semi-major  axis,  a,  and  eccentricity,  e. 

For  the  purpose  of  the  discussion  in  the  next  section  on  the  PAR,  we  will  consider  equality  constraints  on  the  semi¬ 
major  axis  and  eccentricity.  Given  a  fixed  a,  £  in  Eq.  (5)  can  be  replaced  with 


(6) 


resulting  in 


p2  +  wip  +  F(p)  +  ^  =  0. 


(7) 


Energy,  is  related  to  eccentricity  (and  orbital  angular  momentum)  by  [21] 


2 


(8) 


However,  note  that  the  angular  momentum  h  is  related  to  the  parameters  (a,  a ,  8 , 8.  p.  p).  One  can  show  that  h  can  be 
expressed  as 


h  =  hip  +  h2p2  +  hop  +  h4 


(9) 


where  (note  that  our  expressions  here  differ  from  those  of  DeMars  and  Jah  [12]) 


hi  =  q  x  Up 

h2  =  Up  x  (d  cos  <5ua  +  (5u^) 

h3  =  Up  x  q  +  q  x  (a  cos  <5ua  +  ju^) 

h4  =  q  x  q. 


Therefore,  the  magnitude  squared  of  the  angular  momentum  is  given  by 


||h||2  =  cop2  +  P{p)p  +  U(p), 


(10) 


where 


P(p)  =  Cip2  +  c2p  +  c3 
U  ( p )  =  C4p4  +  C5p3  +  Cop 2  +  c7p  +  Cg, 
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and  where  Cj  (i  =  0, . . . ,  8)  are  given  by 

|i,  m2 

co  =  ||hi|| 
ci  =  2hi  ■  I12 
c2  =  2hi  ■  h3 
c3  =  2hi  ■  h4 

|i,  ii2 

C4  =  ||h2|| 

C5  =  2h2  •  h3 

c6  =  2h2  •  h4  +  ||h3||2 

C7  =  2h3  •  h4 

ii,  1 1 2 

cs  =  ||h4j| 

Substituting  Eq.  (10)  into  Eq.  (8),  which  in  turn  gets  substituted  into  Eq.  (5),  one  obtains 

04/5^  +  03j6^  +  02p2  +  Gtl/3  +  Clo  =  0,  (11) 

where  a,  (1  =  0. ....  4)  are  functions  of  p  and  e  and  are  given  by: 

ao  =  F(p)U{p)  +  p2{  1  -  e2) 
ai  =  F(p)P(p)  +  wiU(p) 
a  2  =  U(p)  +  CqF(p)  +  wiP(p) 
a3  =  P(p)  +  c0wi 
a4  =  Cq. 


For  given  values  of  a  and  e,  in  addition  to  the  attributable  variables  (a,  d,  <5, 5)  one  can  then  solve  the  two  nonlinear 
equations  Eq.  (7)  and  Eq.  (11)  for  the  unknowns  p  and  p.  It  is  doubtful,  however,  that  there  exists  a  closed  form 
solution  to  these  equations. 


3.  PROBABILISTIC  ADMISSIBLE  REGION 

The  general  goal  is  to  obtain  the  probabilistic  characterization  of  the  uncertainty  in  the  variables  (p.  p)  given  knowledge 
of  the  statistics  of  the  angles-only  measurement  process.  In  other  words,  we  would  like  to  obtain  a  probability  density 
function  (pdf)  p(p,  p)  that  characterizes  the  PAR. We  do  so  as  follows.  Firstly,  notice  that  one  can  view  the  least  squares 
solution  for  (a,  d,  5,  8)  from  the  short-arc,  angles-only  measurements  as  a  mapping  through  which  one  can  map  the 
uncertainty  in  measurements  to  uncertainty  in  (a,  d,  5,  8)  at  epoch.  Let  that  pdf  be  denoted  by  p(a ,  d,  <5, 8).  Further 
assume  that  we  have  a  probabilistic  assessment  of  the  distribution  over  the  semi-major  axis,  a,  and  the  eccentricity,  e, 
with  pdfs  p(a)  and  p(e),  respectively.  We  assume  that  the  semi-major  axis  and  the  eccentricity  are  independent  of  each 
other  and  both  from  (a,  d,  <5,  i5)*.  Distributions  in  a  and  e  can  be  obtained  from  any  known  information  and  physics 
of  the  RSO  population.  In  other  words,  the  joint  distribution  in  (a,  d,  8,  6.  a,  e )  is  given  by: 

p(a,  d,  S,  S,  a,  e)  =  p(a,  d,  S,  S)p(a)p(e). 

Using  Eq.  (7)  and  Eq.  (11)  one  can  map  p(a,  d,  6 , 5,  a,  e)  to  obtain  p(p,  p). 

Such  an  uncertainty  map  is  impossible  to  achieve  in  closed-form,  even  when  the  measurement  process  is  Gaussian. 
This  is  due  to  the  nonlinearity  present  in  Eq.  (7)  and  Eq.  (11).  Therefore,  in  the  present  paper  we  will  use  the  Monte 
Carlo  method  to  obtain  a  particle  representation  of  the  uncertainty  in  (p,p).  Note  that  the  mapping  of  Gaussian 
measurements  to  angles  and  angle-rates  via  the  least  squares  mapping  retains  Gaussianity.  In  other  words  p(a,  d,  S ,  d) 
will  be  Gaussian  under  the  assumption  of  a  Gaussian  measurement  noise  model.  While  in  the  next  section  we  will 
assume  that  the  measurement  process  is  Gaussian,  the  Monte  Carlo  procedure  described  in  the  next  section  still  holds 
as  long  as  we  can  sample  the  pdf  that  describes  the  measurement  process.  If  we  are  not  able  to  sample  the  distribution, 
approximate  MC  methods,  e.g.,  importance  sampling,  can  be  employed  instead. 

*In  general,  this  is  not  true.  The  measurement  process  statistics  are  dependent  on  the  orbital  geometry,  e.g.,  a  and  e. 
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4.  A  MONTE-CARLO  REPRESENTATION  OF  THE  PAR 


Assuming  that  we  know  the  observation  process  noise  characteristics,  Based  on  the  discussion  presented  in  the  pre¬ 
vious  section,  one  can  in  principle  map  the  uncertainty  in  the  angles-only  measurement  process  to  the  uncertainty  in 
their  rates  given  the 

It  is  assumed  that  the  sensor  noise  characteristics  are  known.  For  the  sake  of  simplicity,  we  assume  that  the  measure¬ 
ment  noise  is  Gaussian  with  means  and  covariances  for  right  ascension  and  declination,  respectively, 

at  time  t j,  i  =  0, . . . ,  n  —  1.  Following  a  batch  least  squares  method  [22],  the  measurements  (a,,  Si),  i  =  0, . . . ,  n  —  1, 
and  their  statistics  can  be  mapped  to  a  Gaussian  distribution  in  (ao,  do ,  S o,  do ) .  Clearly,  we  can  sample  this  distribu¬ 
tion  for  the  MC  analysis.  Let  the  particles  in  these  variables  be  denoted  by  {(a^\  cxq\  S{f\  S^)},  j  =  1 ,N, 
where  N  is  the  chosen  number  of  particles  to  represent  the  uncertainty  cloud.  Note  that  by  our  assumption  that  the 
semi-major  axis  and  eccentricity  are  independent  of  the  angles  and  the  angle-rates,  sampling  of  the  the  latter  vari¬ 
ables  can  be  done  by  directly  sampling  from  p(cto,  do,  So,  So).  Likewise,  sample  p(a)  directly  to  obtain  semi-major 
axis  particles  {a^)},  and  p(e)  to  obtain  eccentricity  particles  (e^)},  j  =  1, . . . ,  N.  Notice  that  due  to  the  statistical 
independence  property,  the  ordering  of  the  particles  inside  each  of  the  three  sets  {(o^,  cxq\  Sq\  Sq'1)},  {aM^ }  and 
}  is  irrelevant,  but  each  particle  from  each  set  has  to  be  uniquely  and  consistently  matched  with  one  from  each  of 
the  other  two  sets  to  obtain  joint  particles  {jd-d  j  =  [a^'1 ,  d[ p ,  Sjf^ ,  Sq^  ,  a*dj .  e(i)  yj  =  1; . . . ;  j\T. 

Next,  each  particle  yd)  can  now  be  inserted  in  Eq.  (7)  and  Eq.  (1 1)  to  obtain  a  particle  =  (p^\  p(J) )  in  range, 
range-rate  space.  The  resulting  set  of  particles  {Z^  }  =  {(p^\  p^)}  is  then  the  desired  particle  cloud  that  represents 
the  PAR. 


5.  CONVERTING  THE  PAR  MC  PARTICLE  CLOUD  INTO  A  PAR  GMM 


While  an  MC  particle  cloud  is  a  faithful  representation  of  the  true  uncertainty  (within  the  bounds  of  the  assumptions 
made  in  the  procedure  described  above),  it  is  difficult  to  perform  many  analytical  computations  (e.g.,  filtering  using 
GMM  techniques)  with  it  in  closed-form.  This  motivates  the  need  to  convert  the  PAR  particle  cloud  into  a  PAR  GMM, 
to  name  one  important  analytical  model.  This  has  been  done  previously  by  DeMars  and  Jah  [12],  where  the  CAR  was 
modeled  using  a  uniform  distribution  and  a  GMM  approximation  of  the  uniform  distribution  was  proposed.  Now  that 
we  have  a  non-uniform  PAR,  one  needs  to  resort  to  more  general  methods  for  representing  the  PAR  using  a  GMM. 
The  standard  procedure  in  the  probabilistic  inference  literature  is  the  EM  method  that  can  convert  a  sample  (i.e.,  a 
particle  cloud)  into  a  GMM  [19].  We  will  not  describe  the  details  of  the  method  here.  Instead  we  refer  the  reader 
to  the  algorithm  in  Ref.  [20],  which  is  a  robust  and  versatile  algorithm  called  the  FJ-EM.  The  method  allows  for  the 
user  to  specify  a  maximum  number  of  GMM  components  and  it  selects  the  number  of  components  that  best  represent 
the  particle  cloud.  The  user  is  also  not  required  to  choose  a  good  initial  guess  of  the  components.  In  other  words, 
the  method  does  not  require  a  careful  initialization  of  the  algorithm  (other  EM  algorithms  do  require  a  very  careful 
initialization.) 


As  an  illustration,  consider  the  2000  particle  cloud  generated  from  the  four  component  planar  GMM  with  weights: 


means: 


and  covariances: 


uh  =  0.25,  i  =  1,2, 3, 4, 


Ml  = 

[10.0 

0.0]T 

m2  = 

[-10. 

0  o.of 

M  3  = 

[0.0  10.0]T 

Mi  = 

[0.0 

-  lO.of 

2.0 

0.0  ' 

•  _ 

0.0 

2.0 

1  ^  — 
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The  generated  particle  cloud  was  then  fed  into  the  FJ-EM  algorithm  and  the  following  5-component  GMM  initial 
guess  was  used  for  initializing  the  algorithm: 


wl  =  0.1 
w°2  =  0.4 
w%  =  0.2 
wl  =  0.2 

w%  =  0.1 


Hi  =  [1.0  0.0]T 
=  [-1.0  0.0]T 
Hi  =  [0.0  i.o] T 

H°4  =  [-1.0  -  1.0]T 
Mb  =  [-5-0  —  5.0]t 


The  algorithm  resulted  in  the  following  set  of  GMM  weights: 


w{  =  0.25 
w[  =  0.24849 
w{  =  0.26156 
w{  =  0.23995 
w[  =  0.0 


where  we  notice  that  one  component  has  been  eliminated.  The  final  set  of  means  were  found  to  be 

n{  =  [9.96428  0.02116]T 
n{  =  [-9.87450  -  0. 02471] T 
h{  =  [0.03449  10.06634] T 
n{  =  [0.04656  -  9. 96678] T 

and  the  final  set  of  covariances  were  found  to  be 

f  _  [  1.92951  0.02156  ' 

^  ~  [  0.02156  1.92172 

/  _  T  2.10782  -0.04989  ' 

“  [  -0.04989  1.92587 

/  _  [  1.89439  0.05116  ' 

^3  “  [  0.05116  2.35834 

/  _  [  1.99178  0.01623  ' 

“  0.01623  1.85653 


The  resulting  GMM  is  graphically  shown  in  Figure  1  against  the  particles.  Note  that  a  component  has  a  weight  of  zero 
(i.e.,  eliminated)  resulting  in  an  effective  number  of  components  of  4.  While  the  true  values  of  the  weights,  means  and 
covariances  were  not  recovered  exactly,  the  error  in  these  parameters  is  quite  small.  In  the  next  section,  we  implement 
the  FJ-EM  algorithm  in  the  PAR  analysis,  converting  the  particle  cloud  in  the  p  —  p  space  into  a  GMM  representation 
of  the  PAR. 
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Fig.  1.  An  example  showing  the  original  means  (equally  weighted)  that  generated  the  shown 
particle  cloud.  The  cloud  was  then  fed  into  the  FJ-EM  algorithm  that  generated  the  GMM 
approximation  shown.  The  FJ-EM  algorithm  was  able  to  very  accurately  reconstruct  the  true 
GMM  that  generated  the  particle  cloud. 


6.  SIMULATION  RESULTS 

In  this  section  we  provide  a  numerical  example  that  demonstrates  the  above  ideas.  We  consider  an  orbit  with  the  orbital 
element  values  shown  in  Table  1  at  the  initial  simulation  time.  Using  the  Socorro,  NM,  ground  sensor,  9  right  ascension 
and  declination  observations  were  collected  at  the  rate  of  one  observation  every  20  seconds.  The  measurement  noise 
is  assumed  to  be  Gaussian  with  an  angular  standard  deviation  of  2  arcsec  for  both  right  ascension  and  declination. 


Table  1.  Orbital  Elements  of  the  True  Orbit 


Orbital  Element 

Value 

Semimajor  Axis  (km) 

26  571. 

Eccentricity 

0.2 

Inclination  (deg) 

55.0 

Perigee  (deg) 

-120.0 

Right  Ascension  of  the  Ascending  Node  (deg) 

-13.24 

Initial  True  Anomaly  (deg) 

110.0 

The  9  measurements  and  their  statistics  are  then  used  to  obtain  the  (Gaussian)  statistics  for  the  attributable  vector  at 
epoch  using  the  method  of  Maruskin  et  al.  [7].  Because  each  of  the  measurement  statistics  is  assumed  Gaussian  and 
the  batch  least-squares  solution  preserves  Gaussianity,  the  attributable  vector  was  found  to  have  the  following  mean: 


A ^oc 

—0.32184 

A t'di 

3.37094 

VS 

1.33430e-4 

.  V/j  . 

8.08795e-5 
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and  covariance: 


6.21086e-11 
—  1.45300e-12 
0.0 
0.0 


—  1.45300e-12 
5.27577e"14 
0.0 
0.0 


0.0 

0.0 

6.21086e-11 
—  1.45300e-12 


0.0 

0.0 

— 1.45300e-12 
5.27577e-14 


The  reason  the  covariance  matrix  has  a  block  form  is  that  the  right  ascension  and  declination  angular  measurements  are 
independent.  Semi-major  axis  is  assumed  be  uniformly  distributed  between  10,000  km  and  50,000  km.  Eccentricity 
is  assumed  to  be  uniformly  distributed  between  0.0  and  0.4. 


For  each  particle  y^\  Eq.  (7)  can  be  used  to  constmct  a  curve  in  the  p  —  p  space  corresponding  to  the  sampled 
value.  Repeating  for  all  particles  y<:,)  we  obtain  a  “cloud”  of  semi-major  axis  curves  (which  we  shall  call  the 
“a-curves”).  Figure  2  shows  the  cloud  generated  from  a  set  of  500  particles  y<J> .  Likewise,  Eq.  (11)  can  be  used  to 
construct  a  cloud  of  curves  (which  we  shall  call  the  “e-curves”),  shown  in  Figure  3.  Recall  that  there  is  exactly 
one  a-curve  that  corresponds  to  one  e-curve.  The  intersection  (if  an  intersection  exists)  of  these  two  sets  of  curves 
results  in  the  particle  cloud  in  the  p  —  p  space  (we  do  not  plot  the  particles,  but  show  the  curves  and  their  intersections 
in  Figure  4). 
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Fig.  2.  The  “cloud”  of  a-curves. 
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Fig.  3.  The  “cloud”  of  e-curves. 


We  selected  a  relatively  small  number  of  particles  above  to  demonstrate  the  non-uniformity  of  the  a-  and  e-curves 
(otherwise,  the  plot  is  too  dense  to  be  able  to  see  the  distribution  of  the  curves).  Using  a  sample  of  4000  particles 
instead  of  500  we  get  the  p  —  p  particle  cloud  sown  in  Figure  5.  The  figure  also  shows  the  GMM  approximation  of 
the  particle  cloud  generated  from  the  FJ-EM  algorithm.  For  the  initial  GMM  we  grid  the  p  —  p  space  into  a  set  of  10 
ten  grid  points  along  each  direction,  for  a  total  of  100  grid  points).  This  grid  defines  the  means  of  the  initial  GMM. 
All  components  had  the  same  weight  of  1/100  =  0.01  and  all  have  the  same  covariances  with  a  standard  deviation  of 
one-tenth  of  an  earth  radius  in  the  p  direction  and  200  meter  per  second  in  the  p  direction.  The  generated  GMM  had 
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Fig.  4.  The  “cloud”  of  a-  and  e-curves  overlaid  on  top  of  each  other. 


a  total  of  38  non-zero  components  (using  a  threshold  of  0.0001  for  the  weights).  As  can  be  seen,  the  GMM  seems  to 
reflect  the  cloud  rather  faithfully.  This  GMM  can  now  be  used  along  the  lines  suggested  by  DeMars  and  Jah  [12]  in  a 
Bayesian  filtering  scheme. 
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Fig.  5.  The  particle  “cloud”  in  p  —  p  space  at  its  FJ-EM-based  GMM  approximation. 


7.  CONCLUSION 

In  this  paper  we  used  a  Monte  Carlo  approach  to  show  how  to  construct  a  probabilistic  admissible  region  in  p  —  p 
space.  The  PAR,  as  shown,  is  clearly  non-uniform.  We  also  proposed  to  use  the  FJ-EM  algorithm  to  convert  the 
particle  cloud  into  a  GMM  representation  of  the  PAR.  This  GMM  can  then  be  used  to  initialize  a  Bayesian  filter.  In 
this  paper  we  addressed  the  situation  where  only  two  constraints  exist  (semi-major  axis  and  eccentricity).  Current 
work  focuses  on  allowing  additional  constraints,  such  as  a  minimum  constraint  on  perigee. 
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